
C    Copyright 2011 Prof K.Sridharan
C    This file is part of SAP2
C
C    SAP2 is free software: you can redistribute it and/or modify
C    it under the terms of the GNU General Public License as published by
C    the Free Software Foundation, either version 3 of the License, or
c    (at your option) any later version.
C
C    SAP2 is distributed in the hope that it will be useful,
C    but WITHOUT ANY WARRANTY; without even the implied warranty of
C    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C    GNU General Public License for more details.
C
C   You should have received a copy of the GNU General Public License
C    along with SAP2.  If not, see <http://www.gnu.org/licenses/>.


	SUBROUTINE PMPYES(QPMP,HPMP,EFF,HSHUT,NPCH,QQ,HH,ETA,WH,WB,FKNR)



	DIMENSION XI(100),WWH(100),WWB(100),
     *            WHS(89,3),WBS(89,3),RMS(100)

	!ms$if .not. defined(LINKDIRECT)
	!ms$attributes dllexport :: PMPYES
	!ms$endif


	real*4 QQ(*),HH(*),ETA(*),WH(*),WB(*),HSHUT,FKNR
	integer NPCH


	PI = 3.14159

	open (unit = 16, file = 'yes1.res')
	open (unit = 17, file = 'yes2.res')
	open (unit = 18, file = 'yes3.res')

        DATA WHS/

     *  0.634,0.643,0.646,0.640,0.629,0.613,0.595,0.575,0.552,0.533,
     *  0.516,0.505,0.504,0.510,0.512,0.522,0.539,0.559,0.580,0.601,
     *  0.630,0.662,0.692,0.722,0.753,0.782,0.808,0.832,0.857,0.879,
     *  0.904,0.930,0.959,0.996,1.027,1.060,1.090,1.124,1.165,1.204,
     *  1.238,1.258,1.271,1.282,1.288,1.281,1.260,1.225,1.172,1.107,
     *  1.031,0.942,0.842,0.733,0.617,0.500,0.368,0.240,0.125,0.011,
     *  -.102,-.168,-.255,-.342,-.423,-.494,-.556,-.620,-.655,-.670,
     *  -.670,-.660,-.655,-.640,-.600,-.570,-.520,-.470,-.430,-.360,
     *  -.275,-.160,-.040,0.130,0.295,0.430,0.550,0.620,0.634,

     * -.690,  -.599,  -.512,  -.418,  -.304,  -.181,  -.078,  -.011,
     *   .032,   .074,   .130,   .190,   .265,   .363,   .461,   .553,
     *    .674,   .848,  1.075,  1.337,  1.629,  1.929,  2.180,  2.334,
     *   2.518,  2.726,  2.863,  2.948,  3.026,  3.015,  2.927,  2.873,
     *   2.771,  2.640,  2.497,  2.441,  2.378,  2.336,  2.288,  2.209,
     *   2.162,  2.140,  2.109,  2.054,  1.970,  1.860,  1.785,  1.571,
     *   1.357,  1.157,  1.016,   .927,   .846,   .744,   .640,   .500,
     *    .374,   .191,   .001,  -.190,  -.384,  -.585,  -.786,  -.972,
     *  -1.185, -1.372, -1.500, -1.940, -2.160, -2.290, -2.350, -2.350,
     *  -2.230, -2.200, -2.130, -2.050, -1.970, -1.895, -1.810, -1.730,
     *  -1.600, -1.420, -1.130,  -.950,  -.930,  -.950, -1.000,  -.920,
     *   -.690,

     *	-2.230,-2.000,-1.662,-1.314,-1.089,-0.914,-0.750,-0.601,-0.440,
     *	-0.284,-0.130, 0.055, 0.222, 0.357, 0.493, 0.616, 0.675, 0.680,
     *	 0.691, 0.752, 0.825, 0.930, 1.080, 1.236, 1.389, 1.548, 1.727,
     *	 1.919, 2.066, 2.252, 2.490, 2.727, 3.002, 3.225, 3.355, 3.475,
     *	 3.562, 3.604, 3.582, 3.540, 3.477, 3.321, 3.148, 2.962, 2.750,
     *	 2.542, 2.354, 2.149, 1.909, 1.702, 1.506, 1.310, 1.131, 0.947,
     *	 0.737, 0.500, 0.279, 0.082,-0.112,-0.300,-0.505,-0.672,-0.797,
     *	-0.872,-0.920,-0.949,-0.960,-1.080,-1.300,-1.500,-1.700,-1.890,
     *	-2.080,-2.270,-2.470,-2.650,-2.810,-2.950,-3.040,-3.100,-3.150,
     *	-3.170,-3.170,-3.130,-3.070,-2.960,-2.820,-2.590,-2.230/


	DATA WBS/ 	

     *  -.648,-.547,-.414,-.292,-.187,-.105,-.053,0.012,0.042,0.097,
     *  0.156,0.227,0.300,0.371,0.444,0.522,0.596,0.672,0.738,0.763,
     *  0.797,0.837,0.865,0.883,0.886,0.877,0.859,0.838,0.804,0.758,
     *  0.703,0.645,0.583,0.520,0.454,0.408,0.370,0.343,0.331,0.329,
     *  0.338,0.354,0.372,0.405,0.450,0.486,0.520,0.552,0.579,0.603,
     *  0.616,0.617,0.606,0.582,0.546,0.500,0.432,0.360,0.288,0.214,
     *  0.123,0.037,-.053,-.161,-.248,-.314,-.372,-.580,-.740,-.880,
     * -1.00,-1.12,-1.25,-1.37,-1.49,-1.59,-1.66,-1.69,-1.77,-1.65,
     *  -1.59,-1.52,-1.42,-1.32,-1.23,-1.10,-.980,-.820,-.684,

     *  -1.420, -1.328, -1.211, -1.056,  -.870,  -.677,  -.573,  -.518,
     *   -.380,  -.232,  -.160,   .000,   .118,   .308,   .442,   .574,
     *    .739,   .929,  1.147,  1.370,  1.599,  1.839,  2.080,  2.300,
     *   2.480,  2.630,  2.724,  2.687,  2.715,  2.688,  2.555,  2.434,
     *   2.288,  2.110,  1.948,  1.825,  1.732,  1.644,  1.576,  1.533,
     *   1.522,  1.519,  1.523,  1.523,  1.490,  1.386,  1.223,  1.048,
     *    .909,   .814,   .766,   .734,   .678,   .624,   .570,   .500,
     *    .407,   .278,   .146,   .023,  -.175,  -.379,  -.585,  -.778,
     *  -1.008, -1.277, -1.560, -2.070, -2.480, -2.700, -2.770, -2.800,
     *  -2.800, -2.760, -2.710, -2.640, -2.540, -2.440, -2.340, -2.240,
     *  -2.120, -2.000, -1.940, -1.900, -1.900, -1.850, -1.750, -1.630,
     *  -1.420,

     *	-2.260,-2.061,-1.772,-1.465,-1.253,-1.088,-0.921,-0.789,-0.632,
     *	-0.457,-0.300,-0.075, 0.052, 0.234, 0.425, 0.558, 0.630, 0.621,
     *	 0.546, 0.525, 0.488, 0.512, 0.660, 0.850, 1.014, 1.162, 1.334,
     *	 1.512, 1.683, 1.886, 2.105, 2.325, 2.580, 2.770, 2.886, 2.959,
     *	 2.979, 2.962, 2.877, 2.713, 2.556, 2.403, 2.237, 2.080, 1.950,
     *	 1.826, 1.681, 1.503, 1.301, 1.115, 0.960, 0.840, 0.750, 0.677,
     *	 0.604, 0.500, 0.352, 0.161,-0.040,-0.225,-0.403,-0.545,-0.610,
     *	-0.662,-0.699,-0.719,-0.730,-0.810,-1.070,-1.360,-1.640,-1.880,
     *	-2.080,-2.270,-2.470,-2.650,-2.810,-2.950,-3.040,-3.100,-3.150,
     *	-3.170,-3.200,-3.160,-3.090,-2.990,-2.860,-2.660,-2.260/

         QQR=QPMP
         HHR=HPMP
         ETAR=EFF	 

	 WRITE (16,3150)
         XI(1) = 45.0
         WWH(1) = HSHUT/HHR
         WH(45) =WWH(1)

         DO I = 2,NPCH
         VDASH=QQ(I)/QQR
         HDASH=HH(I)/HHR
         EDASH=ETAR/ETA(I)
         BETA=VDASH*HDASH *EDASH
         WWH(I)=HDASH/(1+VDASH**2)
         WWB(I)=BETA/(1+VDASH**2)
         DD=PI/ATAN(VDASH)
         XI(I)=44/DD+45
         WRITE (16,3100)QQ(I),HH(I),ETA(I),VDASH,HDASH,EDASH,BETA,
     *   WWH(I),WWB(I),DD,XI(I)
         END DO
3100     FORMAT (2X, F6.3, 2X, F6.1, 2X, F5.1, 6(1X,F7.4),2(1X,F6.2))
3150     FORMAT (3X,'Q',6X,'H',6X,'ETA',5X,'v=Q/QR',2X,'h=H/HR',3X,
     *   'nr/n',4X,'Beta',4X,'WH',7X,'WB',5X,'d',5X,'Ix')

         IJ = 46
         YY=IJ
	 K=1
         DO WHILE (YY.LE.XI(NPCH))

         DO I =1,NPCH-1
         IF (YY.GT.XI(I).AND.YY.LE.XI(I+1)) K=I
         END DO
	 K1=K-1
	 K2=K
	 K3=K+1
         IF(K.EQ.1) THEN
         K1=1
         K2=2
	 K3=3
         ENDIF

	 AA=(WWH(K1)-WWH(K2))/((XI(K1)-XI(K2))*(XI(K1)-XI(K3)))-
     *      (WWH(K2)-WWH(K3))/((XI(K2)-XI(K3))*(XI(K1)-XI(K3)))
	 BB=(WWH(K1)-WWH(K2))/(XI(K1)-XI(K2))-AA*(XI(K1)+XI(K2))
	 CC=WWH(K1)-AA*XI(K1)**2-BB*XI(K1)
         WH(IJ) = AA*YY**2 + BB*YY + CC
         IF (IJ.EQ.56) WH(IJ) = 0.5

         IJ = IJ+1
         YY=IJ
         END DO

         IJ = 45
	 IFLAG=0
	 DO WHILE(IFLAG.EQ.0)
	 IF(XI(2).GT.IJ)THEN
	 IJ=IJ+1
	 ELSE
	 IFLAG=1
	 ENDIF
	 END DO
	 IJST=IJ
         YY=IJ
	 K=1
         DO WHILE (YY.LE.XI(NPCH))

         DO I =1,NPCH-1
         IF (YY.GT.XI(I).AND.YY.LE.XI(I+1)) K=I
         END DO
	 K1=K-1
	 K2=K
	 K3=K+1
         IF(K.EQ.1.OR.K.EQ.2) THEN
         K1=2
	 K2=3
         K3=4
         ENDIF

	 AA=(WWB(K1)-WWB(K2))/((XI(K1)-XI(K2))*(XI(K1)-XI(K3)))-
     *      (WWB(K2)-WWB(K3))/((XI(K2)-XI(K3))*(XI(K1)-XI(K3)))
	 BB=(WWB(K1)-WWB(K2))/(XI(K1)-XI(K2))-AA*(XI(K1)+XI(K2))
	 CC=WWB(K1)-AA*XI(K1)**2-BB*XI(K1)
         WB(IJ) = AA*YY**2 + BB*YY + CC
         IF (IJ.EQ.56) WB(IJ) = 0.5
		             
         IJ = IJ+1
         YY=IJ
         END DO

         IJEND = IJ -1
	 DO II=1,3
	 RMS(II)=0
	 END DO
         DO I = 45, IJEND
	 DO II=1,3
         RMS(II)= RMS(II)+(WH(I) - WHS(I,II))**2
         END DO
	 END DO
	 NN=IJEND-45+1
	 RMSMIN=1000
	 DO II=1,3
	 RMS(II)=SQRT(RMS(II)/NN)
	 IF(RMS(II).LT.RMSMIN)THEN
	 ISPCH=II
	 RMSMIN=RMS(II)
	 ENDIF
	 END DO

	 IF(ISPCH.EQ.1)THEN
	  FKNR = 0.7
	 ELSEIF(ISPCH.EQ.2)THEN
	  FKNR = 2.2
	 ELSE
	  FKNR = 1.1
	 END IF


         DIFF = WH(45) - WHS(45,ISPCH)
         DO I = 35,44
         WH(I) = WHS(I,ISPCH) + DIFF*0.1*(I-35)
         END DO

	 IGAP=IJST-35
	 DIFF=WB(IJST)-WBS(IJST,ISPCH)
	 DO I=35,IJST-1
	 WB(I)=WBS(I,ISPCH)+DIFF*(1.0/IGAP)*(I-35)
	 END DO

         DIFF = WH(IJEND) - WHS(IJEND,ISPCH)
         DIFF1 = WB(IJEND) - WBS(IJEND,ISPCH)
         DO I = IJEND+1, IJEND + 5
         WH(I) = WHS(I,ISPCH) + DIFF*0.2*(IJEND+5-I)
         WB(I) = WBS(I,ISPCH) + DIFF1*0.2*(IJEND+5-I)
         END DO

         DO I = 1,34
         WH(I) = WHS(I,ISPCH)
         WB(I) = WBS(I,ISPCH)
         END DO
         DO I = IJEND+6,89
         WH(I) = WHS(I,ISPCH)
         WB(I) = WBS(I,ISPCH)
         END DO
	
c	ENDIF

	DO I=1,89
	WRITE(17,3300)I,WH(I),WHS(I,ISPCH)
	END DO
3300	FORMAT(2X,I5,2F10.3)
	DO I=1,89
	WRITE(18,3300)I,WB(I),WBS(I,ISPCH)
	END DO

	close(16)
	close(17)
	close(18)

	end




	SUBROUTINE PMPNO(ISPCH,WH,WB)


	!ms$if .not. defined(LINKDIRECT)
	!ms$attributes dllexport :: PMPNO
	!ms$endif


	real*4 WH(*),WB(*)
	integer ISPCH


	DIMENSION WHS(89,3),WBS(89,3)



	open (unit = 16, file = 'no1.res')
	open (unit = 17, file = 'no2.res')
	open (unit = 18, file = 'no3.res')








        DATA WHS/

     *  0.634,0.643,0.646,0.640,0.629,0.613,0.595,0.575,0.552,0.533,
     *  0.516,0.505,0.504,0.510,0.512,0.522,0.539,0.559,0.580,0.601,
     *  0.630,0.662,0.692,0.722,0.753,0.782,0.808,0.832,0.857,0.879,
     *  0.904,0.930,0.959,0.996,1.027,1.060,1.090,1.124,1.165,1.204,
     *  1.238,1.258,1.271,1.282,1.288,1.281,1.260,1.225,1.172,1.107,
     *  1.031,0.942,0.842,0.733,0.617,0.500,0.368,0.240,0.125,0.011,
     *  -.102,-.168,-.255,-.342,-.423,-.494,-.556,-.620,-.655,-.670,
     *  -.670,-.660,-.655,-.640,-.600,-.570,-.520,-.470,-.430,-.360,
     *  -.275,-.160,-.040,0.130,0.295,0.430,0.550,0.620,0.634,

     * -.690,  -.599,  -.512,  -.418,  -.304,  -.181,  -.078,  -.011,
     *   .032,   .074,   .130,   .190,   .265,   .363,   .461,   .553,
     *    .674,   .848,  1.075,  1.337,  1.629,  1.929,  2.180,  2.334,
     *   2.518,  2.726,  2.863,  2.948,  3.026,  3.015,  2.927,  2.873,
     *   2.771,  2.640,  2.497,  2.441,  2.378,  2.336,  2.288,  2.209,
     *   2.162,  2.140,  2.109,  2.054,  1.970,  1.860,  1.785,  1.571,
     *   1.357,  1.157,  1.016,   .927,   .846,   .744,   .640,   .500,
     *    .374,   .191,   .001,  -.190,  -.384,  -.585,  -.786,  -.972,
     *  -1.185, -1.372, -1.500, -1.940, -2.160, -2.290, -2.350, -2.350,
     *  -2.230, -2.200, -2.130, -2.050, -1.970, -1.895, -1.810, -1.730,
     *  -1.600, -1.420, -1.130,  -.950,  -.930,  -.950, -1.000,  -.920,
     *   -.690,

     *	-2.230,-2.000,-1.662,-1.314,-1.089,-0.914,-0.750,-0.601,-0.440,
     *	-0.284,-0.130, 0.055, 0.222, 0.357, 0.493, 0.616, 0.675, 0.680,
     *	 0.691, 0.752, 0.825, 0.930, 1.080, 1.236, 1.389, 1.548, 1.727,
     *	 1.919, 2.066, 2.252, 2.490, 2.727, 3.002, 3.225, 3.355, 3.475,
     *	 3.562, 3.604, 3.582, 3.540, 3.477, 3.321, 3.148, 2.962, 2.750,
     *	 2.542, 2.354, 2.149, 1.909, 1.702, 1.506, 1.310, 1.131, 0.947,
     *	 0.737, 0.500, 0.279, 0.082,-0.112,-0.300,-0.505,-0.672,-0.797,
     *	-0.872,-0.920,-0.949,-0.960,-1.080,-1.300,-1.500,-1.700,-1.890,
     *	-2.080,-2.270,-2.470,-2.650,-2.810,-2.950,-3.040,-3.100,-3.150,
     *	-3.170,-3.170,-3.130,-3.070,-2.960,-2.820,-2.590,-2.230/


	DATA WBS/ 	

     *  -.648,-.547,-.414,-.292,-.187,-.105,-.053,0.012,0.042,0.097,
     *  0.156,0.227,0.300,0.371,0.444,0.522,0.596,0.672,0.738,0.763,
     *  0.797,0.837,0.865,0.883,0.886,0.877,0.859,0.838,0.804,0.758,
     *  0.703,0.645,0.583,0.520,0.454,0.408,0.370,0.343,0.331,0.329,
     *  0.338,0.354,0.372,0.405,0.450,0.486,0.520,0.552,0.579,0.603,
     *  0.616,0.617,0.606,0.582,0.546,0.500,0.432,0.360,0.288,0.214,
     *  0.123,0.037,-.053,-.161,-.248,-.314,-.372,-.580,-.740,-.880,
     * -1.00,-1.12,-1.25,-1.37,-1.49,-1.59,-1.66,-1.69,-1.77,-1.65,
     *  -1.59,-1.52,-1.42,-1.32,-1.23,-1.10,-.980,-.820,-.684,

     *  -1.420, -1.328, -1.211, -1.056,  -.870,  -.677,  -.573,  -.518,
     *   -.380,  -.232,  -.160,   .000,   .118,   .308,   .442,   .574,
     *    .739,   .929,  1.147,  1.370,  1.599,  1.839,  2.080,  2.300,
     *   2.480,  2.630,  2.724,  2.687,  2.715,  2.688,  2.555,  2.434,
     *   2.288,  2.110,  1.948,  1.825,  1.732,  1.644,  1.576,  1.533,
     *   1.522,  1.519,  1.523,  1.523,  1.490,  1.386,  1.223,  1.048,
     *    .909,   .814,   .766,   .734,   .678,   .624,   .570,   .500,
     *    .407,   .278,   .146,   .023,  -.175,  -.379,  -.585,  -.778,
     *  -1.008, -1.277, -1.560, -2.070, -2.480, -2.700, -2.770, -2.800,
     *  -2.800, -2.760, -2.710, -2.640, -2.540, -2.440, -2.340, -2.240,
     *  -2.120, -2.000, -1.940, -1.900, -1.900, -1.850, -1.750, -1.630,
     *  -1.420,

     *	-2.260,-2.061,-1.772,-1.465,-1.253,-1.088,-0.921,-0.789,-0.632,
     *	-0.457,-0.300,-0.075, 0.052, 0.234, 0.425, 0.558, 0.630, 0.621,
     *	 0.546, 0.525, 0.488, 0.512, 0.660, 0.850, 1.014, 1.162, 1.334,
     *	 1.512, 1.683, 1.886, 2.105, 2.325, 2.580, 2.770, 2.886, 2.959,
     *	 2.979, 2.962, 2.877, 2.713, 2.556, 2.403, 2.237, 2.080, 1.950,
     *	 1.826, 1.681, 1.503, 1.301, 1.115, 0.960, 0.840, 0.750, 0.677,
     *	 0.604, 0.500, 0.352, 0.161,-0.040,-0.225,-0.403,-0.545,-0.610,
     *	-0.662,-0.699,-0.719,-0.730,-0.810,-1.070,-1.360,-1.640,-1.880,
     *	-2.080,-2.270,-2.470,-2.650,-2.810,-2.950,-3.040,-3.100,-3.150,
     *	-3.170,-3.200,-3.160,-3.090,-2.990,-2.860,-2.660,-2.260/

        

        DO I=1,89
        WH(I)=WHS(I,ISPCH)
        WB(I)=WBS(I,ISPCH)
        END DO

	DO I=1,89
	WRITE(17,3300)I,WH(I),WHS(I,ISPCH)
	END DO
3300	FORMAT(2X,I5,2F10.3)
	DO I=1,89
	WRITE(18,3300)I,WB(I),WBS(I,ISPCH)
	END DO

	close(16)
	close(17)
	close(18)


	END
